function R2 = IPCA_oos(R,Z,IsF,X,W,Nt,K,varargin)

% R is N*T
% Z is N*T*L
% IsF is N*T
% X is L*T
% W is L*L*T
% PSF is M*T

% Standard IPCA call as IPCA_oos(R,Z,IsF,X,W,Nt,K)
% IPCA with PSF call as IPCA_oos(R,Z,IsF,X,W,Nt,K,PSF)
% IPCA with ONLY PSF call as IPCA_oos(R,Z,IsF,X,W,Nt,0,PSF)

if nargin>7 && ~isempty(varargin)
    PSF = varargin{1};
    PSF_version = true;
else
    PSF_version = false;
end

[N, T] = size(R);
L = size(X,1);

THalf = round(T/2);

numer1_i = NaN(N,T); numer2_i = NaN(N,T); denom_i = NaN(N,T);
numer1_x = NaN(L,T); numer2_x = NaN(L,T); denom_x = NaN(L,T);
Tangency_F = NaN(T,1); Tangency_PSF = NaN(T,1);
F_F = NaN(T,K);

for t = THalf:T-1
    if PSF_version
        [Gamma, F] = IPCA(X(:,1:t),W(:,:,1:t),Nt(1:t),K,PSF(:,1:t));
    elseif ~PSF_version
        [Gamma, F] = IPCA(X(:,1:t),W(:,:,1:t),Nt(1:t),K);
    end
    
    Zt = squeeze(Z(:,t,:));
    Wt = W(:,:,t);
    Xt1 = X(:,t+1);
    Rt1 = R(:,t+1);
    ValidIndex = IsF(:,t+1);

    Exp1_i = zeros(N,1); Exp2_i = zeros(N,1);
    Exp1_x = zeros(L,1); Exp2_x = zeros(L,1);

    if K>0
        GammaBeta = Gamma(:,1:K);
        Ft1 = (Gamma'*Wt*Gamma)\(Gamma'*Xt1);
        LambdaF = mean(F,2,'omitnan');
        SigmaF = cov(F','omitrows');
        Betat = Zt*GammaBeta;
        Exp1_i = Exp1_i + Betat*Ft1;
        Exp2_i = Exp2_i + Betat*LambdaF;
        Exp1_x = Exp1_x + Wt*GammaBeta*Ft1;
        Exp2_x = Exp2_x + Wt*GammaBeta*LambdaF;
        
        w = SigmaF\LambdaF; w = w/sum(w);
        Tangency_F(t+1) = w'*Ft1;
        F_F(t+1,:) = Ft1;
    end

    if PSF_version
        GammaDelta = Gamma(:,K+1:end);
        PSFt1 = PSF(:,t+1);
        LambdaPSF = mean(PSF(:,1:t),2,'omitnan');
        SigmaPSF = cov(PSF(:,1:t)','omitrows');
        Deltat = Zt*GammaDelta;
        Exp1_i = Exp1_i + Deltat*PSFt1;
        Exp2_i = Exp2_i + Deltat*LambdaPSF;
        Exp1_x = Exp1_x + Wt*GammaDelta*PSFt1;
        Exp2_x = Exp2_x + Wt*GammaDelta*LambdaPSF;
        
        w = SigmaPSF\LambdaPSF; w = w/sum(w);
        Tangency_PSF(t+1) = w'*PSFt1;
    end
    
    Err1t1 = Rt1 - Exp1_i;
    Err2t1 = Rt1 - Exp2_i;
    numer1_i(ValidIndex,t+1) = Err1t1(ValidIndex);
    numer2_i(ValidIndex,t+1) = Err2t1(ValidIndex);
    denom_i(ValidIndex,t+1) = Rt1(ValidIndex);

    Err1t1 = Xt1 - Exp1_x;
    Err2t1 = Xt1 - Exp2_x;
    numer1_x(:,t+1) = Err1t1;
    numer2_x(:,t+1) = Err2t1;
    denom_x(:,t+1) = Xt1;
end
%numer1_x(end,:) = []; numer2_x(end,:) = []; denom_x(end,:) = []; % the last 'managed' portfolio is just the average return; remove it

Alpha = mean(numer1_i,2,'omitnan'); Rbar = mean(denom_i,2,'omitnan'); RelPrcErr_i = sum(Alpha.^2,'omitnan')/sum(Rbar.^2,'omitnan');
numer1_i = numer1_i.^2; numer2_i = numer2_i.^2; denom_i = denom_i.^2;
TotalR2_i = 1 - sum(numer1_i(:),'omitnan')/sum(denom_i(:),'omitnan');
PredR2_i = 1 - sum(numer2_i(:),'omitnan')/sum(denom_i(:),'omitnan');

R2i = 1 - sum(numer1_i,2,'omitnan')./sum(denom_i,2,'omitnan'); 
Ti = sum(isfinite(numer1_i),2); 

TimeSeriesR2_i = sum(Ti.*R2i,'omitnan')/sum(Ti,'omitnan');


R2x = 1 - sum(numer1_i,1,'omitnan')./sum(denom_i,1,'omitnan'); CrossSectionR2_i = mean(R2x,'omitnan');

Alpha = mean(numer1_x,2,'omitnan'); Rbar = mean(denom_x,2,'omitnan'); RelPrcErr_x = sum(Alpha.^2,'omitnan')/sum(Rbar.^2,'omitnan');
numer1_x = numer1_x.^2; numer2_x = numer2_x.^2; denom_x = denom_x.^2;
TotalR2_x = 1 - sum(numer1_x(:),'omitnan')/sum(denom_x(:),'omitnan');
PredR2_x = 1 - sum(numer2_x(:),'omitnan')/sum(denom_x(:),'omitnan');
R2i = 1 - sum(numer1_x,2,'omitnan')./sum(denom_x,2,'omitnan'); Ti = sum(isfinite(numer1_x),2); TimeSeriesR2_x = sum(Ti.*R2i,'omitnan')/sum(Ti,'omitnan');
R2x = 1 - sum(numer1_x,1,'omitnan')./sum(denom_x,1,'omitnan'); CrossSectionR2_x = mean(R2x,'omitnan');

%--------------------------------------------------------------------------
R2.TotalR2_i = TotalR2_i;
R2.PredR2_i = PredR2_i;
R2.TotalR2_x = TotalR2_x;
R2.PredR2_x = PredR2_x;
R2.TimeSeriesR2_i = TimeSeriesR2_i; 
R2.CrossSectionR2_i = CrossSectionR2_i;
R2.TimeSeriesR2_x = TimeSeriesR2_x;
R2.CrossSectionR2_x = CrossSectionR2_x;
R2.RelPrcErr_i = RelPrcErr_i;
R2.RelPrcErr_x = RelPrcErr_x;
R2.F_F = F_F;

if K>0
    R2.TangencySR_F = sqrt(12)*mean(Tangency_F,'omitnan')/std(Tangency_F,'omitnan');
end
if PSF_version
    if mean(Tangency_PSF,2,'omitnan') < 0
        Tangency_PSF = -Tangency_PSF;
    end
    R2.TangencySR_PSF = sqrt(12)*mean(Tangency_PSF,'omitnan')/std(Tangency_PSF,'omitnan');
end
return